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ABSTRACT 

We report new observations of the far infrared and submillimctcr water vapor emission of NML 
Cygnus based on data gathered with the Infrared Space Observatory and the Submillimeter Wave 
Astronomy Satellite. We compare the emission from NML Cyg to that previously published for VY 
CMa and W Hya in an attempt to establish the validity of recently proposed models for the outflow 
from evolved stars. The data obtained support the contention by Ivezic & Elitzur (1995, 1997) that 
the atmospheres of evolved stars obey a set of scaling laws in which the optical depth of the outflow is 
the single most significant scaling parameter, affecting both the radiative transfer and the dynamics 
of the outflow. Specifically, we provide observations comparing the water vapor emission from NML 
Cyg, VY CMa and W Hya, and find, to the extent permitted by the quality of our data, that the 
results are in reasonable agreement with a model developed by Zubko & Elitzur (2000). Using this 
model we derive a mass loss based on the dust opacities, spectral line fluxes, and outflow velocities 
of water vapor observed in the atmospheres of these oxygen-rich giants. For VY CMa and NML Cyg 
we also obtain an estimate of the stellar mass. 

Subject headings: circumstellar matter - infrared: stars - stars: mass loss - stars: winds, outflows - 
stars: individual (NML Cygnus, VY Canis Majoris, W Hydrae) 



1. INTRODUCTION 

NML Cyg, also designated IRC+10448, has been studied extensively since its discovery four decades ago in the near 
infrared sky survey conducted by Neugebauer, Martz & Leighton (1965). Though the star is quite faint at visual 
wavelengths, it is enormously luminous in the infrared. Morris and Jura (1983) associated this star with the Cyg OB2 
association, which placed the star at ~ 2 kpc and defined its luminosity as 5 x 10 5 £©. These authors concluded 
that the star has a mass of 50 M and a minimum mass loss of ~ 6 x 10~ 5 Mq yr~ x , indicated by the dense dust 
shell surrounding the star. The star exhibits significant variability in the infrared. Monnier et al (1997) report flux 
variations of ±0.25 mag in the N band at 10.2 /im with an irregular period around 940 days. 

Using speckle interferometry in the near infrared in the K, L, and M bands, Ridgway et al. (1986) estimated the 
optical depth at 3.4/im to be 2.4 and derived an inner radius of the dust shell of ~ 0.045 arcsec, where they estimated 
the dust temperature to be > 1200 K. Through an extensive de-reddening analysis these authors also found consistency 
with the previously proposed spectral type M6 III (Wing, Spinrad & Kuhi, 1967) and concluded that the star's effective 
temperature is T e // ~ 3250 K. More recently, Monnier et al. (1997) conducted long baseline interferometry at 11/im 
and modeled the dust outflow in somewhat greater detail, obtaining better fits to all gathered data by assuming a 
stellar effective temperature of only T e ff ~ 2500 K. 

The star's line of sight velocity relative to the local standard of rest is approximately known. Cohen, et al (1987) 
investigated the 1612 MHz OH-maser emission from NML Cyg. They found the distribution among the many individual 
maser features they were able to resolve at 0.06 km s _1 to be double peaked with peaks centered roughly at a velocity 
with respect to the local standard of rest vi sr = ±22 km s _1 , suggesting that the central star's line of sight velocity 
is close to zero. These authors assumed the star to lie at a distance of 2 kpc, which placed the radius of the maser 
emitting shell around 7.3 x 10 16 cm from the star. The star's radius, based on its color temperature and luminosity at 
that distance, was then ~ 250 times smaller, ~3x 10 14 cm. Diamond, Norris and Booth (1984), who extended these 
OH observations, confirmed a similar stellar velocity. 
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A similar study concentrating on the water vapor masers associated with outflow from NML Cyg by Richards et al. 
(1996) shows a rather different picture, exhibiting a strong component at vi sr = —20 km s -1 , with no other emission 
besides a hint of a second feature at ~ +6 km s -1 . This asymmetry is explained by the authors as due to a collimated 
outflow close to the star at an angular distance of ~ 100 mas. This outflow, they conjecture, gradually becomes less 
collimated at a larger distance of order a few arcseconds, where the OH 1612 MHz masers are observed. Observations 
carried out at several epochs indicate proper motion of order 20 km s _1 , assuming a distance of 2 kpc. This would 
be consistent both with the radial velocity of the OH maser peaks and the approaching velocity of the H 2 maser 
emission. 

A paper by Danchi et al. (2001) reports on the dust distribution interferometrically measured at 11/im in the 
infrared. They conclude that there are two distinct dust shells, both of which have moved away from the star by 
approximately the same distance over a period of ~ 6 yr. Their separation suggests that they were formed ~ 65 yr 
apart. If the transverse outflow is comparable to the radial velocities of the masers, these authors favor a distance to 
the star of ~ 1200 pc, rather than 2 kpc, preferred by Richards et al. (1996). 

The dust emission measured by Danchi et al. (2001) peaks at roughly 100 mas from the star but is distributed 
rather broadly over distances from ~ 25 to > 150mas. Determining a displacement of only ~ 25 mas of the diffuse 
dust outflow, measured over the span of only 6 yr, appears to us to be a more difficult observation than measuring 
the angular displacement of the more compact H2O masers over a period of ~ 13 yr, with quite small error bars. For 
purposes of our current paper, we therefore assume the distance of 2 kpc for NML Cyg, suggested by Richards et al. 
(1996). 

Boboltz & Marvel (2000) have examined the v = 1,J = 1 — ► SiO maser emission around the star with the Very 
Long Baseline Array. They find that the masers form a roughly elliptical ring around the star, with dimensions 
~ 8 x 10 14 by 10 15 cm. The maser velocity distribution together with the spatial distribution led these authors to 
conclude that the SiO maser shell rotates with a velocity wsini ~ 11 km s , about a systemic velocity of —6.6 km 
s _1 . While it is difficult to see how such a high rotational velocity could be induced at such a large distance from 
the star, a comparable system velocity emerges from consideration only of the double-peaked distribution of maser 
velocities, with peaks at ~ 5 and —18 km s _1 , averaging out to ~ —6.5 km s _1 . 

The star's velocity can also be derived from its thermal SiO v = emission. This has a roughly parabolic spectral 
form, suggesting spherically symmetric outflow. Lucas et al, (1992) found that the peak emission occurs at about -5 
km s _1 , and the mean emission velocity is about -1 km s _1 . These authors also determined the radius at the half 
power point for SiO emission to lie at ~ 2 x 10 16 cm, well beyond the SiO maser emitting region but also well inside 
the shell of OH masers. They derived a mass loss of order 1.6 x 10~ 4 M Q for the star, somewhat higher, but not at 
significant variance with the estimates of Morris and Jura (1983). 

Many of these properties of NML Cyg arc summarized in Table 1 , where comparative values for two other oxygen-rich 
stars, VY Canis Majoris and W Hydrae, can also be found. 



Table 1. Characteristics of VY CMa, W Hydrac fc NML Cyg 



Star 


VY CMa a 


W Hydrac 6 


NML Cyg 


Star's Distance (pc) 


1500 


115 


2000*" 


Star's Luminosity (Lq) 


5 x 10 5 


11,050 


5 x 10 5 


Star's Temperature (K) 


2800 


2500 


2500 d 


Star's Radius (cm) 


2.25 x 10 14 


3.9 x 10 13 


2.6 x 10 14 d 


Final Outflow Velocity (km s _1 ) 


20 


10 


~ 25 c 


Optical Depth 


A., ~ 3.2" 


0.83 at 5500A 


~ 2 at ll^im 6 


Gas-to-dust ratio 


~ 100 


850 


~ 100 6 


Mass loss rate Mq yr -1 


~ 2 x 10~ 4 


2.3 x 10~ 6 


~ 2 X 10~ 4 b 



a Data from Neufeld et al. (1996) and Barlow et al. (1996) 
b For sources of data see text 
c NML Cyg data, this paper 



d Model of Monnier et al. (1997) 
e optical depth in the J band 



Since NML Cyg is one of the most luminous sources in the northern sky it was extensively used as a calibration source 
for the Short Wavelength Spectrometer (SWS) Fabry-Perot (F-P) mode, during the Infrared Space Observatory (ISO) 
mission. Some of the spectral data reported here were obtained as part of these calibration procedures. On ISO the star 
was also studied as part of several dedicated programs aimed at understanding the atmospheres of evolved oxygen-rich 
stars. One of these programs concentrated on detecting the water vapor emission from such stars with the aid of the 
Submillimeter Wave Astronomy Satellite (SWAS) as well as with SWS and the Long Wavelength Spectrometer (LWS) 
on ISO. ISO data were, respectively, obtained for W Hya by Barlow et al. (1996) and Neufeld et al. (1996), and for 
VY CMa by Neufeld et al. (1996, 1999), while SWAS data for these two stars were published by Harwit & Bergin 
(2002). Taken together, these data tended to support the theoretical models for the atmospheres and outflow from 
oxygen rich stars developed by Ivezic and Elitzur (1995, 1997) and Zubko & Elitzur (2000). 

Section 2 of the present paper reports water vapor observations obtained with ISO. Section 3 reports similar obser- 
vations obtained with SWAS. Section 4 briefly describes the theoretical models for stellar outflow devised by Zubko & 
Elitzur (2000). Section 5 makes use of their model to obtain a mass-loss rate for NML Cyg. Section 6 compares the 
observed water-vapor line strengths to those that the model predicts. Section 7 discusses some apparent discrepancies 



3 



in the data and attempts to plausibly account for them. In section 8 we discuss the balance between gravitational 
attraction and repulsion by light pressure, which may account for some of the apparent anomalies in outflow velocities 
and mass-loss rates. 

2. H 2 OBSERVATIONS WITH ISO 

Early water vapor observations of NML Cyg were obtained by Justtanont et al. (1996) who used the SWS in grating 
mode and, in addition to other spectral features, detected rovibrational transitions of H2O at 2.7 and 6.2 fim and 9 
pure rotational lines of the molecule between 36 to 45 /im. The resolving power of the grating was of order 680 at 
30^m — insufficient to clearly determine line shapes. From continuum observations and an assumed outflow velocity 
of 27.7 km s _1 these authors estimated a dust mass loss of ~ 2 x 1CP 6 M Q yr _1 , which translates into an overall mass 
loss of ~ 2 x 1CP 4 Mq yr _1 for a dust-to-gas mass ratio of 100. 

The ISO observations cited in the present paper were carried out with the telescope pointed at 20 ft 46 m 25.46 s , +40° 
6' 59.6" (J 2000). With the SWS/F-P, water vapor features were obtained at a spectral resolving power of ~ 30, 000. 
Most of these mid-infrared lines emanate from high-lying energy levels of H 2 and consistently exhibit P Cygni profiles, 
clearly seen at the higher spectral resolution of the Fabry-Perot instrument. The line shapes are consistent with the 
assumption that the lines are emitted close to the star where the stellar outflow is still at a high temperature, and 
that cooler H2O, moving at considerably higher outflow velocities at greater distances from the star, absorbs much of 
the emitted light. Figure 1 shows several of the spectral lines registered with the SWS/F-P, and shows the line fits 
that we adopted. 

At a resolving power of 30,000 the SWS/F-P emission lines are well resolved. In order to recover the emission line 
strength before self-absorption by cool water vapor in the outer parts of the outflow, we assumed that the source 
spectrum consists of a continuum flux F c plus a parabolic emission feature F — Fq ■ (v — vq) 2 /v 2 for the range 
\v — vq\ < w e , where vq is the line of sight velocity of the star and v e is the expansion velocity in the emitting region. 
This velocity is expected to be different for different emission lines, and provides information on the expansion velocity 
in different temperature regimes. 

We then fitted the parabola convolved with the instrument profile to most closely replicate the portion of the emission 
feature on the long-wavelength side of the peak. Since this portion of the emission line should be least affected by self- 
absorption in an accelerating outflow, this procedure permitted us to restore the emitted flux on the short wavelength 
side of the emission peak and estimate the total emitted flux. These steps tended to increase the line flux by 20 to 
30% relative to the directly measured emission component in the P Cygni profiles. To fully simulate the P Cygni 
profile, we further postulated a constant fractional absorption / a (, s over some velocity range v m i n to v max applied to 
the continuum plus fully restored parabolic emission. All three parameters, f a bs, ^min and v max were allowed to be 
free. For a given value of f a b s we selected v m i n to give the best fit to the short-wavelength side of the absorption 
feature and chose v max to give the best fit to the short-wavelength side of the emission feature and the long- wavelength 
side of the absorption feature. Between these three parameters we generally obtained a good fit. The peaks of the 
emission lines tend to fall around vi sr = —10. We find v max to be of order —50 km s _1 for most of our fits, while the 
star's line of sight velocity with respect to the local standard of rest is taken to be vi sr ~ km s _1 . Such a high value 
for v max is somewhat surprising, as discussed in Section 7, below. 

The P Cygni profiles in the mid-infrared water vapor features of NML Cyg are consistent with similar profiles 
observed in VY CMa (Neufeld et al., 1999), while the same lines observed by Neufeld et al. (1996) in W Hya show no 
self-absorption, indicating a significantly lower column depth in the star's cool outer layers. 

The spectral lines observed by the LWS/F-P were obtained at a spectral resolving power of ~ 10,000. These lines 
involve lower-lying energy states emitted by cooler gas at greater distances from the star, and exhibit no noticeable 
self absorption. Presumably this is due to their emission in the more tenuous outer regions of the outflow where 
self- absorption is less likely. Figure 2, shows four emission lines obtained with the LWS/F-P. The data were reduced 
with the ISO LWS/F-P interactive data reduction program. At 95 /xm the resolving power of the instrument is ~ 9800 
and gradually declines to about 9500 at 125 /im. 

To obtain these data we repeatedly stepped the LWS Fabry-Perot instrument through a wavelength range that 
spanned roughly 7 to 8 steps on cither side of the line which generally encompassed another 4 to 6 steps. These roughly 
20 to 22 steps were equally spaced. Each step was displaced approximately half a spectral resolution element or ~ 15 
km s _1 from the previous step. The entire spectral range was stepped through repeatedly in what was designated 
as the "rapid-scanning mode" of operation. In reducing the data, a sloping spectral base line was subtracted. The 
reduced spectra are shown in Figure 2. The best-fit spectral peak is displaced to shorter wavelengths relative to the 
laboratory values, with a consistent velocity displacement of 12 — 14 km s _1 , in good agreement with the SWS/F-P 
values. Full widths at half maximum range from ~ 31 to ~ 69 km s _1 . Because the spectral resolving power in this 
wavelength range, ~ 30 km s _1 , is comparable to the observed line widths, and significantly broader than the spectral 
velocity resolution of the SWS spectra or the SWAS spectra reported in Section 3, we fitted the observed emission 
features with a single Gaussian shape plus a baseline, a choice primarily determined by the Fabry-Perot transmission 
function. This line fit exhibits broader wings than a parabolic fit might suggest, but appears more appropriate under 
the circumstances and, in any case, does not significantly affect the results we report. 

3. H 2 OBSERVATIONS WITH SWAS 

NML Cyg was observed in the 557 GHz li — loi (538^m) ground state ortho transition with SWAS. The spectral 
resolution with this instrument is of order 3 x 10 5 , i.e. about 1 km s _1 . The observed line is considerably weaker 
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than the infrared lines, but nevertheless was detectable. The data were gathered in observations covering 165 hours 
during two different observing sessions. The first covered 106 hours spreading over the period from November 8, 2002 
to December 20, 2002. The second observing session covered 59 hours and took place between April 26 and May 15, 
2003. The emission line together with a parabolic fit is shown in Figure 3. The spectral resolving power in this spectral 
range is more than an order of magnitude higher than in the SWS and LWS ranges, and the observed line shape is 
dominated entirely by the actual outflow. The parabolic fit we have chosen reflects the expectation that the outflow 
may, to first order, be isotropic in the star's rest frame. The flux integrated over the line is 1.3 x 10~ 20 W cm -2 , the 
peak flux is 184 Jy and the peak of the emission line lies at vi sr = km s^ 1 . The one sigma uncertainty of the fitted 
integrated flux is 0.18 x 10~ 20 W cm" 2 . 

The emission observed in this spectral range is due to cold water vapor. The upper level of this transition lies only 
T = hvjk ~ 26. 7K above the ground ortho-H 2 state, so that it is readily excited in collisions with cool gas in the 
outer regions of the star's outflow. 

4. THEORETICAL MODELS OF THE OUTFLOW FROM EVOLVED STARS 

Until the past few years, models of the outflow from evolved stars generally treated the radiative transfer through 
the outflowing atmosphere separately from the dynamics of the flow. This tended to lead to conflicting notions: On the 
radiative transfer side, the outflow was considered to proceed at some constant terminal velocity. From the dynamic 
perspective, dust grains were pictured as absorbing starlight which they quickly re-emitted in the infrared while also 
transferring absorbed momentum to the gas, thus accelerating the outflow. But if the infrared radiation is so closely 
tied to momentum transfer, it cannot simultaneously be made consistent with a constant outflow velocity — leading 
to the cited conflict. 

This difficulty was first overcome in a series of papers by Ivezic & Elitzur (1995, 1997) and Zubko & Elitzur (2000), 
who produced a self-consistent outflow model combining both dust radiative transfer and hydrodynamics. In oxygen- 
rich evolved stars, water vapor is a dominant coolant in the outflow. The model exhibits scaling laws that define the 
radial variations in the velocity and the radiation fields in terms of a radial distance measured in units equal to the 
inner diameter of the dust absorption shell. The gas in the outflow is largely heated through collisions with the grains 
that also transfer linear momentum to the gas, thereby accelerating it. Over a significant portion of the outflow the 
balance between the heating and cooling is largely dominated by water vapor. The H2O is mainly heated by absorption 
of radiation emanating from dust in the immediate vicinity of the star, though some heating also is due to collisions 
with hydrogen molecules that have been collisionally heated by dust grains. H 2 cannot rapidly cool through radiative 
transitions since it is restricted to quadrupole emission at the temperatures found in outflows. Most of the heating 
and cooling takes place through pure rotational H2O transitions throughout the bulk of the outflow, the cooling being 
almost directly proportional to the heating throughout the flow (Zubko & Elitzur, private communication, 2001). 
As discussed by Harwit & Bergin (2002) this nearly constant proportionality has the effect of simulating a radiative 
transfer and temperature profile of a gray body, scaling the temperature T to radial distance R from the star through 
the constancy of the product T 4 R 2 . 

5. DERIVED MASS LOSS FROM NML CYG AND VY CMA 

Figure 4 shows representative points from the infrared continuum spectrum of NML Cyg obtained through ground- 
based photometry (Strecker & Ney 1974), with additional data from IRAS. These points are fitted with a synthetic 
spectrum modeled on the procedure of Zubko & Elitzur (2000), who assumed a silicate- rich dust envelope in the 
outflow from oxygen-rich stars and used the DUSTY code of Ivezic, Nenkova, & Elitzur (1999). Figure 5 provides 
matching data for VY CMa. The observational data are obtained from Harwit et al. (2001). Input parameters for the 
computations that produce the fitted spectra of Figures 4 and 5 are shown for the respective stars in Table 2. 

The spectral fits to the observed dust continuum spectra shown in Figures 4 and 5, while not ideal, seem satisfactory, 
particularly given the high variability of each star. The temperature and outflow velocity profile derived with the model 
of Zubko & Elitzur is shown for NML Cyg and VY CMa in Figure 6. Again, it is worth noting the roughly constant 
value of Tg as i? 2 , already mentioned in Section 4. 

The outflow models constructed here for NML Cyg and VY CMa are essentially different from the one established 
by Zubko & Elitzur (2000) for W Hydrae, in that the wind is optically thick for both NML Cyg and VY CMa, with 
t v ~ 70 to 80, whereas for W Hydrae, t v ~ 1. This means that the dust envelope defines the total flux in most of 
the spectrum, as indicated in Figures 4 and 5 by the merger, beyond ~ 10/zm, of the dotted line representing dust 
emission and the solid line fitted to the total observed flux. 

Table 2 underlines the similarities between NML Cyg and VY CMa. One major difference, however, is the ten times 
higher H 2 abundance in VY CMa, required by the star's large observed H 2 line fluxes. 

To check the stability of our results for NML Cyg, we varied the temperature of the central star over the range 
2,500 ± 400K, the terminal outflow velocity over the range 25 ± 10km s _1 , and the optical depth t v over values of 
75 ± 10, constrained largely by the near infrared portion of the spectrum in the 1-8 fxm region. Combining these 
variations, with dust temperatures 1400 ± 100 K at the inner edge of the outflow, well within the realm constrained by 
the far infrared spectrum in the 20 to 100 /iin region, and a ratio of outer to inner edge of the outflow Y = 1, 000 ± 100, 
we derived a range of outflow rates of M — 7 ± 2 x 10~ 4 M Q yr _1 and a gas-to-dust ratio 350 ± 30. Our estimates 
of M and final velocity, however, are subject to an additional uncertainty of ~ 30%, and the gas-to-dust ratio to an 
additional uncertainty of ~ 60%, due to the way that the DUSTY model takes gravity into account. This is further 
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Tabic 2. Input, Fitting, and Derived Parameters for NML Cyg and VY CMa 

NML Cyg VY CMa 



Input Parameters: 

Stellar BB temperature T* 

Luminosity L- M 

Distance 

Terminal outflow velocity a,fa 
Fitting Parameters: 

Optical depth r v at 0.55/^m 
Dust temperature at i?i nner 

Y = Houter /dinner 

H2O abundance at -Ri nne r 
Ortho-H 2 : para-H 2 
Derived Parameters: 
Mass loss rate M 
Gas-to-dust mass ratio 
Stellar radius 
Shell inner radius -Ri llncr 



2500 K 
5 x 1O 5 L 

2000 pc 
25 km s _1 

75 
1400 K 
1000 

4 x 10~ 5 
1 



6.9 x 10 

350 
2.6 x 10" 
1.05 x 10 1 



2800 K 
5 x 10 5 L o 

1500 pc 
20 km s _1 

75 
1400 K 
1000 
4 x 10~ 4 
1 

9.0 x 10~ 4 
510 

2.1 x 10 14 cm 
1.15 x 10 15 cm 



Silicate dust optical data from Laor & Draine (1993) 

Dust size distribution from Mathis, Rumple & Nordsicck (1977) 



discussed in section 8. Given the quality of available data, these results appear to be satisfactorily robust. The 
uncertainties in the fitting and derived parameters for VY CMa are approximately the same. 

Values cited in the "Model" columns of Table 3 provide our best predictions of integrated line fluxes, especially for 
spectral lines not yet observed. Future observations thus will be able to pass judgement on the general validity of the 
model. 



Table 3. Relative H 2 Q Flux Levels from W Hydrac, VY CMa and NML Cyg 



Wavelength 


Transition 




T( 


W Hydrae 


VY CMa 




NML Cyg 


A (Aim) 










(K) 


(K) 


10~ 20 Wcm~ 2 


10~ 19 Wcm~ 2 


10~ 20 Wcm~ 2 










Observed 


Model" 


Observed" 


Model d 


Observed 


Model" 


538.289 


o: 


lio 




loi 


61.0 


34.2 


0.45 


0.76 


0.34 


0.57 


1.3 


1.2 


179.527 


o: 


2l2 




loi 


114 


34.2 


8.66 


12.5 


8.7 


5.7 




16.5 


174.624 


o: 


3o3 




2l2 


196 


114 


9.21 


8.3 


8.8 


3.8 




10.7 


125.356 


P 


4(14 




3l3 


319 


204 


16.6 


12.4 




5.2 


14.6 


14.1 


113.538 


o: 


4l4 




3o3 


323 


196 


17.4 


11.3 




5.6 


10.1 


14.7 


108.073 


o: 


2l2 




lio 


194 


61.0 


13.0 


15.5 


13.2 


8.5 


17.5(7?) 


19.1 


100.983 


p: 


220 




111 


196 


53.4 




22.7 




10.6 


32.9 


21.8 


99.493 


o: 


505 




4l4 


468 


323 




12.6 




5.6 


40.0 


15.3 


95.627 


p: 


5l5 




4 4 


469 


319 




12.2 




5.2 


14.3 


15.1 


89.989 


P 


322 




2n 


297 


137 


27.3 


18.1 




8.8 


24.5(7?) 


20.1 


78.742 


o: 


4 2 3 




3l2 


432 


249 


28.0 


14.1 




7.4 


37.8(7?) 


17.7 


66.438 


o: 


330 




221 


411 


194 


22.9 


18.5 


17.2(7?) 


10.3 




21.6 


58.699 


o: 


432 




321 


550 


305 


19.5 


15.4 




8.9 


52.5(7?) 


18.8 


40.760 


P 


633 




524 


952 


599 




22.4 


11.7(7?) 


11.3 


53 


24.9 


e 40.691 


o: 


432 




303 


550 


196 


23.0 


36.1 


36.9(7?) 


20.2 




43.6 


37.566 


P 


744 




633 


1335 


952 




13.7 




8.3 


37 


15.8 


36.212 


P 


624 




5l5 


867 


469 




33.7 


5.6(7?) 


17.2 


48 


36.7 


35.938 


o: 


652 




541 


1279 


878 




13.9 


12.1(7?) 


9.0 


47 


18.3 


"31.772 


o: 


441 




3l2 


702 


249 


63.0 


40.0 


24.1 


24.6 




47.5 


29.837 


o: 


725 




616 


1126 


643 


32.0 


30.7 


23.5 


19.8 




39.6 



a ISO from Neufcld ct al. (1996), Barlow ct al. (1996); SWAS data, this paper and Harwit & Bcrgin (2002) 
b W Hydrac model from Zubko & Elitzur (private communication, 2001) 
c For sources of data see text 
d This paper 

e Blending may need to be taken into account 

(??) Question marks indicate less reliably established line fluxes 



6. COMPARISON OF THEORY AND OBSERVATIONS 

Once a best-fit outflow model was in hand, we again used the Zubko & Elitzur (2000) model to derive the expected 
line strengths for all the water vapor lines, based on a best fit to the relatively few actually observed water vapor 
features, respectively, in NML Cyg, VY CMa and W Hydrae. 

Because the model does not take into account the self absorption of emitted radiation in the stellar outflow, the 
SWS emission lines exhibiting P Cygni profiles have to be adjusted, as described in Section 2, to recover the originally 
emitted flux level. In Table 3 we present the water vapor emission line data for NML Cyg, W Hya, and VY CMa 
obtained with ISO and SWAS, along with the model reductions that Zubko & Elitzur developed for W Hya (Harwit & 



6 



Bergin, 2002) and similar new calculations for NML Cyg and VY CMa. The model assumes that scaling laws apply, 
and we do find rough self-consistency in that respect. Table 3 shows that the water vapor spectra of the three stars, 
despite their considerable differences in luminosity and distance appear to approximately scale relative to each other, 
VY CMa providing line fluxes that are roughly ten times higher than W Hya and about four times higher than NML 
Cyg. 

7. SOME APPARENT ANOMALIES 

Somewhat surprising in the data on NML Cyg is the rather large displacement of the absorption trough seen in 
the P Cygni profiles of mid-infrared lines originating in regions of high temperature surrounding the star. As Figure 
1 shows, some of the inferred outflow velocities in the absorbing gas reach values of the order of 40 km s _1 , though 
with a mean value more nearly 30 km s _1 . These velocities are higher than the displacement of the far-infrared lines 
observed in the LWS Fabry-Perot mode which uniformly are centered on roughly -15 km s _1 and have full widths at 
half maximum ranging from <~ 30 to 70 km s _1 for different, somewhat noisy line shapes. The SWAS data at 557 GHz 
also suggest lower outflow velocities; the observed line, though weak, appears centered on vi sr ~ with a full width 
at half maximum around 40 km s _1 . Comparable velocity profiles were observed by Lucas et al. (1992) for the pure 
rotational J = 2 — > 1 thermal SiO emission at 86.85 GHz. A somewhat similar, though less marked discrepancy in 
line velocities also appears for VY CMa, as judged from the data presented by Neufeld et al. (1999) and Harwit & 
Bergin (2002). 

For NML Cyg, the prime discrepancy in the outflow velocities lies between ISO data, which concentrated on lines 
emitted in higher temperature regions, and the data on lower temperature regions gathered by SWAS and in observa- 
tions of SiO thermal emission. These apparent discrepancies could be resolved if some of the hot absorbing gas near 
the star were streaming out at a substantially higher velocity than cooler gas at larger distances. Diamond, Norris 
and Booth (1984) found similarly puzzling results in their observations of 1612 MHz OH maser emission from NML 
Cyg. They suggested the existence of an inner component with an expansion velocity of 34 km s^ 1 with an extent 
of 1.5 arcsec and an outer component with an expansion velocity of 20 km s _1 at an angular radius of 2.5 arcsec. 
Even the inner of these two regions, however, would lie at a distance of ~ 5 x 10 16 cm from the star and, if in local 
thermal equilibrium, would be far too cool to excite the levels responsible for the P Cygni absorption. Nevertheless, 
it is possible that different expansion components, even in the inner regions surrounding the star, are flowing out 
at different velocities. This could also explain the substantial differences in line widths found for the various LWS 
emission lines. The larger widths could be an artifact of the line-fitting procedure, since the data are quite noisy. 
However, they might also be real, reflecting velocity differences at different radial distances from the star. In future 
searches, it would be useful to see whether substantial velocity differences exist at different radial distances in the 
outflow from NML Cyg. 

Our mass loss rate for NML Cyg and VY CMa is also substantially higher than that cited by other observers. While 
Morris and Jura (1983) only give a minimum mass loss of 6 x 10 _5 M Q yr -1 for NML Cyg, this is an order of magnitude 
lower than the value we give in Table 2. Their calculation is based on the idea of an inverse Stromgren sphere, i.e. the 
ionization of the outer edge of the outflow from NML Cyg by ultraviolet radiation from Cyg OB2 star 5. The authors 
assume this ionizing star to be at the same distance from Earth as NML Cyg, i.e. at the minimum possible distance 
from NML Cyg, ~ 100 pc, given the two stars' angular separation of ~ 2.9°. The radial distance of the ionization 
front from NML Cyg is then assumed to lie on a circumstellar surface at which the ionization rate just balances the 
outflow rate. The mass loss rate computed in this fashion can only be a minimum, and scales as the square of the 
actual distance between the stars. 

Knapp et al. (1982) used observations of the J = 1 - line of CO, by Zuckerman et al. (1977), to determined a 
mass loss of 1.8 x 10~ 4 M Q yr _1 for an outflow velocity of 21 km s _1 , if their assumed stellar distance is changed from 
200 pc to 2 kpc. Their mass loss rate is based on the assumption that the observed CO flux is produced at a surface 
where heating by collisions with ambient gas keeps the CO in thermal equilibrium with that gas. Beyond this radial 
distance and the angular disk that it subtends on the sky, the flux is assumed to drop to zero. The calculated mass 
loss rate then becomes independent of the [CO]/[H2] ratio. 

A third estimate of the mass loss for NML Cyg, due to Netzer & Knapp (1987), is ~ 1.6 x 10~ 4 M Q yr _1 , when 
adjusted to a distance of 2 kpc. It is based on the production of OH in the photodissociation of H 2 by the interstellar 
radiation field at the outer surface of the outflow from the star. The observed OH maser shell radii together with 
the wind outflow velocity then provide a mass loss rate, on the assumption that the abundance of water vapor in the 
outflow is [H20J/H2] = 3 x 10~ 4 . The calculated mass loss also depends on the assumed interstellar radiation field. 
The estimated error bars the authors place on their mass loss rates are plus or minus a factor of 2. However, the 
assumption of a typical interstellar radiation field, appears to be in conflict with the ionizing radiation field assumed 
by Morris & Jura (1983). It is not clear which of the two assumptions, if cither, is correct. 

8. DISCUSSION 

The mass loss of the three stars considered here is generally assumed to be largely driven by the starlight absorbed 
by dust grains. The approach we have taken in this paper is to calculate the mass loss rate derived from the coupled 
dynamic and radiative transfer equations, based largely on dust opacity (Ivezic, Nenkova, & Elitzur, 1999). The 
assumptions we make are thus quite different from those of Moris & Jura (1983), Knapp et al. (1982), or Netzer 
& Knapp (1987). The luminosity L of a star, when totally absorbed by dust, can lead to an absorption of outward 
directed momentum amounting to L/c per second. If the outflow accelerates from a negligibly low initial velocity to 
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high final velocity, v, as indicated in our Figure 6, then the amount of mass that can be accelerated to a velocity v 
in each second is M — L/cv. Taking the assumed luminosity L = 2 x 10 39 erg/sec and a final velocity of 25 km/sec, 
leads to M = 2.67 x 10 22 g s _1 or 4 x 10~ 4 M Q yr _1 . This is of the order the mass loss indicated by the computer 
model we used. 

Our model assumes a spherical outflow and that all of the star's light is absorbed by the dust. That the outflow 
is spherical seems to be in accord with the parabolic line shape that Harwit & Bergin (2002) found for the 557 GHz 
outflow from VY CMa, and the generally parabolic line fits for stars with massive outflows observed for many other 
stars, e.g. in the CO J = 2 - 1 line (Knapp et al., 1982). About the nearly total absorption of star light, there also can 
be little doubt, judging from the optical depths observed for NML Cyg and VY CMa. Two factors, however, do enter 
consideration. The first is the initial outflow velocity of gas just before dust condensation and the onset of radiation 
driven outflow. The second is the gravitational attraction of the star. 

The peak outflow velocity Vi nner of the gas beyond the inner radius where dust is formed is not currently known, 
but it appears to be high. The spectral lines emitted and absorbed in the hottest portions of the outflow, i.e. nearest 
to the star have considerably higher velocities than the flow at larger distances. This can be explained if gravitational 
attraction plays a decisive role once a packet of gas has been accelerated to peak velocity and then is shielded from 
further radiative acceleration by a new dust layer arising between it and the star. 

While the DUSTY computer model of Ivczic, Nenkova, & Elitzur (1999) appears to be the most advanced program 
available at the moment, an inherent limitation is that it takes the force of gravitational attraction into account as 
a reduction in the radiative repulsive force through a parameter T — T va dl Fgrav T is then taken to be constant at 
all distances. This is a valid approximation for optically thin outflows. But when the optical depth of an outflow 
is high, the radiation pressure at large distances from the star drops drastically as dust at large distances from the 
star no longer is able to absorb the far-infrared flux reradiatcd near the inner radius of the outflow. In this case 
the gravitational attraction of the star can become dominant beyond the inner edge of the dusty outflow where the 
radiation pressure drops. We, therefore, may expect a reduction in outflow velocities at large distances from the star, 
a property not reflected in our computer generated outflow velocity in Figure 6. 

If we take the short wavelength outflow velocities Vi nner to be equal to (L/(cM) +i>o, where i>o is the outflow velocity 
at Rinner, the inner edge of the dust cloud where the dust forms, then the final outflow velocity Vf, at large distances 
from the star, determined by the half-width of the 557 GHz line, will be determined by 

L ( 2MG x 1/2 



Vf = V ° + ~Tt ~\ ~R ' (1) 

CM \ riinner J 

For NML Cyg, we can set vj ~ 19 km s _1 , and u + L/(cM) — v inner = 42 km s _1 . This tells us that 2MG / Ri nner <~ 
23km s _1 . Choosing Rinner ~ 10 15 cm then yields the mass of the star as M ~ 20M Q . Cooling of the gas will not 
change this outcome significantly, since the internal energy of gas at T ~ 1400 K, characteristic of the temperature at 
Rinner, is only of order 6 x 10 10 erg g _1 . We can also solve for the mass loss, as 

M ~ (2) 

C{Vinner - V ) 

For NML Cyg, we then have M > 2.4 x 10~ 4 M Q yr _1 , on the assumption that v a must be greater than or equal to 
zero if the gas is to recede sufficiently far from the star to permit dust to form. 

Corresponding velocities for VY CMa are Vf ~ 17km s _1 and Vinner ~ 30km s _1 . As in Table 2 we also take 
Rinner ~ 1-15 x 10 15 cm. These parameters lead to M ~ 8M© and M > 3.3 x 10 4 M Q yr" 1 . 

For both stars we expect v to be significantly lower than v in ner, implying that v in ner is dominated by radiative 
acceleration. This is supported by the findings of Richards et al. (1996) who noted significant acceleration in the H 2 
maser outflows of NML Cyg, suggesting a low value of v . If so, the lower limit to the mass outflow is close to the 
actual value. 

It, therefore, appears that both the total mass loss and the diminishing velocities as a function of distance from the 
star can be roughly accounted for. Within factors of ~ 2, the mass loss rates we derive appear consistent with those 
of other workers, who derived their values from widely differing assumptions. Our mass estimate for VY CMa seems 
low, unless this star has already lost a substantial fraction of its initial mass. Unfortunately, we still are some way 
from narrowing these differences and reaching a satisfactory understanding of massive outflows. 
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Fig. 1. — The spectra and model fits for data from the ISO SWS/F-P instrument. In all panels the wavelength of the H2O transition is 
indicated in the upper left corner. The solid dark lines arc ISO data from the SWS/F-P instrument. The dashed lines represent the two 
components of our model, the absorbed continuum (displaced downward by 1000 Jy in order to be seen more clearly) and the optically 
thick emission from an expanding shell (see text for more details). The dotted lines show the models convolved with instrument response. 
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Fig. 2. — The spectra and model fits for data from ISO LWS. The laboratory wavelengths of the spectral lines are indicated in the box 
on each plot. The dotted lines show the fit with a single Gaussian shape plus a baseline. 
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Fig. 3. — The SWAS data of the ortho H2O lio — ► loi transition at 538.289 fim, smoothed to a channel width of 4.5 km/s. The dotted 
line shows a fit of an optically thick shell with an expansion velocity of 26 km/s. The displayed velocity is with respect to the local standard 
of rest. 




IR photometry • 

IRAS O 

Total emerging flux 

Stellar attenuated flux 

Dust scattering 

Dust emission 



FIG. 4. — The continuum spectrum of NML Cyg. Modeled components are identified in the legend. Photometric data indicated by filled 
circles are from Strecker & Nye (1974). 
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Fig. 5. — The continuum spectrum of VY CMa. Modeled components are identified in the legend. Photometric data indicated by filled 
circles are from Harwit et al. (2001). 
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Fig. 6. — Velocity and temperature profile of the gas flowing out of NML Cyg and VY CMa, derived with the model of Zubko & Elitzur 
(2000). 



